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1.  Introduction 

Three  basic  aspects  of  the  numerical  solution  of  incompressible 
Navier-Stokes  equations  have  been  considered,  namely 

1)  The  question  of  achieving  stable  discretizations 
of  the  incompressibility  constraint. 

2)  The  problem  of  obtaining  accurate  solutions 
in  the  limit  Re  **  go  . 

3)  Devising  of  efficient  numerical  solution 
algorithms  for  solving  the  nonlinear  algebraic 
systems  of  equations  arising  from  the 
discretization  step. 

The  incompressible  Navier-Stokes  equations  are 


du 

dt 


-  vAu  +  (Vu)u  =  - 


div  u  =  0 


~  Vp  +  A  f 
P  P  — 


(1. 


with,  typically,  u  =  0  on  fixed  boundaries  of  the  flow  field  and, 
possibly,  conditions  on  u  specified  at  infinity  (external  flows) . 
In  (1.1),  u  represents  the  velocity  vector,  p  the  pressure, 
p  the  constant  density,  v  the  kinematic  viscosity,  and  f  the 
t>ody  force  per  unit  mass.  For  simplicity,  we  consider  henceforth 

on  <  '>  r'  (I  4  6  Approvsdr^r 
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'  5  r  ■  Informal  tor,  Division 

the  case  of  an  internal  flow,  in  which  the  flow  domain  a  is 
bounded . 

A  standard  weak  form  for  (1.1)  in  the  stationary  case  is 

_il 

as  follows:  find  u  €  Hq (b)  such  that 

V  j*  Vu  :  Vv  -  i  f  u-Vu*v  -  u*Vv*u 

U  ^  a  “ 

■  J  p  d.iv  v  =  f  f-v  V  v  €  H  i  ( il)  (1.2) 

a  U  u 

J  q  div  u  =  0  V  q  €  L^Q) 

a  u 

-.1 

where  Hq(w)  is  t*le  first  order  Sobolev  space  of  vector  functions 

2 

zero  on  the  boundary  and  L^fa)  is  the  class  of  square 

integrable  functions  with  mean  zero  over  U.  The  density  has  been 
normalized  to  p  =  1  in  (1.2). 

Numerical  approximations  to  (u,p)  are  generated  by  the 
following  scheme:  pick  finite  dimensional  subspaces  t*1  ^  Hg(a) 
and  Sh  c:  Lq(h)  and  see^  x  s*1  such  that 

r  n  h  „  h  lphnhh  hnhh 
V  J  vu  :  Vv  -  j  J  u  • Vu  *v  -  u  • Vv  -u 

Q  U 

-  X  phdiv  v*1  =  X  fvh  V  vh  t  H~(Q)  (1.3) 

a  ~  a  u 

X  qhdiv  Uh  =  0  ¥  qh  €  L^(Q) 

a  u 

In  terms  of  this  formulation,  our  work  on  1)  above  deals  with 
stable  choices  of  the  spaces  l/*1  and  S*1  as  summarized  in 


section  2  below.  Work  on  2)  concerns  the  case  v  *•  0  in 
(1.2) -(1.3)  and  is  summarized  in  section  3.  Finally,  the  topic 
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3)  concerns  the  solution  of  the  nonlinear  systems  of  algebriac 
equations  arising  from  (1.3)  and  is  discussed  in  the  fourth 
section. 


2.  Work,  relating  to  discretisation 

This  work  is  a  continuation  of  that  undertaken  during  the 
previous  year's  grant  AFOSR-80-0091 .  Mathematically,  the  following 
condition  is  necessary  for  convergence  of  the  discrete  approxima¬ 
tions  (u\pk)  as  h  -*  0: 

sup  J  qhdiv  vh  >  Yl|qhli  V  qh  fc  Sh  (2.1) 

vhtVh  l* 

|vh|=l 

where  y  >  0  is  independent  of  h.  Previously,  we  gave  a  local 
(elementwise)  test  which  applied  to  all  elements  except  the 
simplest  pairs.  Applying  this  test  reveals  whether  a  given 
element  pair  (U^,S^)  is  stable  or  not.  Unfortunately,  the  simplest 
pairs  of  elements  are  of  considerable  interest  in  practical 
computing,  so  it  is  necessary  to  somehow  prove  or  disprove  that 
(2.1)  is  valid  for  them.  A  major  result  achieved  has  been  to 
show  that  many  often  used  low  order  element  pairs  are,  in  fact, 
unstable  in  the  sense  of  (2.1).  In  addition,  new,  simple  low  order 
element  pairs  which  we  proved  to  be  stable  have  been  introduced. 
Moreover,  a  "postprocessing"  operation  nas  been  given  which  can 
often  be  applied  to  the  unstable  numerical  results  which  restabilize 
them.  This  work  was  carried  out  with  Ph.D.  student  J.  M.  Boland. 

Mr.  (now  Dr.)  Boland  has  been  supported  by  AFOSR  during  his  studies 
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and  graduated  Ph.D.  in  January  1983  from  Carnegie-Mellon  University. 
He  is  presently  Assistant  Professor  of  Mathematics  in  the 
University  of  Pittsburgh.  I  am  confident  he  will  become  an 
outstanding  applied  mathematician. 

The  work  referred  to  in  this  section  is  the  subject  of 
references  [1]  and  [2]  listed  below. 


3.  The  second  topic,  the  inviscid  limit  v  -*  0,  has  lead  to  the 
following  conclusions.  If  h  is  kept  fixed  and  solutions  computed 
for  a  sequence  of  values  of  v  "*  0  then,  as  is  well  known,  the 
numerical  solutions  develop  oscillations  of  nonphysical  character 
of  ever  increasing  amplitude.  Essentially,  this  is  caused  by  the 
fact  that  no  mechanism  exists  in  the  numerical  scheme  for  switching 
off  the  no-slip  boundary  condition,  so  that  for  v  =  0  we  have  an 
overdetermined  system  of  equations.  The  oscillations  are  a  mani¬ 
festation  of  this  potential  overdetermination.  On  the  other  hand, 
if  v  is  kept  fixed,  and  h  -*  0,  then  mathematical  theorems 
are  available  to  quarantee  (subject  to  conditions  of  the  type 
2.1)  convergence  of  numerical  solutions  to  the  true  solution, 

(u,p) .  The  problem  with  letting  h  -*  0  is  the  increasingly 
large  size  of  equation  systems  that  must  be  solved,  if  v  is 
at  all  small.  The  approach  taken  by  the  author  is  based  on  the 
idea  that  h  needs  to  be  small  only  in  certain  locations,  namely 

in  boundary  layers,  and  that,  crucially,  the  thickness  of  such 

1/2 

layers  is  always  0(v  ),  rather  than  0 ( v )  as  would  be  the  case 

for  normal  (as  opposed  to  shear)  layers.  Normal  layers  do  not 
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seem  to  be  an  issue  in  incompressible  flows,  except,  possibly  as 

a  consequence  of  an  artificially  imposed  outflow  condition,  or 

1/2 

some  similar  situation.  The  0(v  )  thickness  means  that 

1/2 

h  =  0(v  )  too  for  stable  solutions  free  of  oscillations, 

and  this  condition  is  not,  in  practice,  too  severe.  This  latter 
approach  is  the  only  one  known  to  the  author  which  actually  solves 
the  posed  problem,  rather  than  some  perturbed  problem,  perturbed 
by  adding  artificial  viscosity  of  amount  far  larger  than  the 
physical  viscosity. 

Work  on  this  topic  continues,  and  will  be  reported  in  [3], 


4.  Solving  the  algebriac  systems  which  arise  from  (1.3)  remains 
a  major  difficulty,  from  the  point  of  view  of  the  amount  of 
computational  labor /cost.  Two  new  approaches  have  been  developed, 
one  dependent  on  time  marching  to  the  steady  state  limit,  and 
the  other  based  on  an  adaptation  of  a  method  used  in  structural 
mechanics  to  the  fluids  case.  To  describe  the  first  method, 
consider  the  iteration 


j~r  (un+1-un)  +  N(un)  +  BPn+1  =  F 


nT  n+1 
B  u  =  g 


n  =  0,1, 


arbitrary. 


(4.1) 


In  (4.1),  the  first  term  is  a  discretization  of  the  second 
denotes  the  viscous  and  convection  terms  evaluated  at  the  previous 
time  level,  and  B  represents  a  discrete  gradient  operator.  The 
fact  that  the  discrete  divergence  operator  is  the  transposed 


•  ‘  t 
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gradient  operator  is  a  property  of  the  continuous  problem  inherited 
by  the  finite  element  formulation. 

T 

To  solve  (4.1),  multiply  the  first  equation  by  B  and  use 

n+1 

the  second  at  n  and  n+1  to  get  the  equation  for  p  , 

BTBpn+1  =  BTF  -  BTN(uk).  (4.2) 

In  (4.2)  the  variable  of  interest  is  Bpn+\  rather  than  pn+^ 
itself,  and  a  special  technique  has  been  devised  for  obtaining 
this  vector,  without  first  finding  pn+^.  Once  this  equation 
is  solved,  un+^  can  be  computed  directly  from  (4.1).  The 
procedure  may  now  be  repeated  to  get  pn+^  and  so  on.  This 
algorithm  has  worked  well  in  practice,  except  that  too  many 
time  steps  are  needed  to  reach  the  steady  state.  This  situation 
is  common  to  most  marching  algorithms. 

The  second,  more  direct  approach  which  has  shown  considerable 
promise  is  based  on  the  following  idea.  The  region  Cl  is  parti¬ 
tioned  into  subdomains  0^,  in  a  way  which  respects  element 
boundaries.  Each  subdomain  will  contain  numerous  elements.  In 
each  subdomain,  the  discrete  Navier-Stokes  equation  is  formed 
by  the  standard  procedure.  The  subdomain  equations  will  interact 
only  across  their  common  boundaries,  and  the  equations  for  a 
given  subdomain  can  be  transformed  in  such  a  way  as  to  express 
the  interior  variables  for  the  subdomain  in  terms  of  its  (unknown) 
boundary  variables.  Linking  the  equations  for  the  boundary 
variables  gives  a  smaller  subsystem  for  solution.  Once  this 
subsystem  is  solved,  the  subdomain  boundary  values  are  known 
and  so  the  subdomain  interior  values  may  be  computed.  This 
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procedure,  which  has  been  used  in  other  fields  for  many  years, 
was  not  used  to  solve  the  Navier-Stokes  equations  because  the 
subdomain  equations  are  singular  (due  to  the  arbitrariness 
in  the  pressure  field)  and  so  the  subdomain  equations  cannot 
£>e  solved.  A  way  of  handling  this  problem  has  been  devised  and 
is  reported  in  [4] .  The  method  is  highly  parallelizable ,  because 
the  subdomain  equations  can  be  processed  simultaneously. 


I  . 
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Papers  prepared  during  grant  period 

1.  Counterexample  to  uniform  stability  of  bilinear-constant, 
velocity-pressure  finite  elements  (with  J.  M.  Boland) 

(to  appear,  Numerische  Mathematik] . 

2.  Stable  and  semistable  low  order  finite  elements  for  viscous 
flows  (with  J.  M.  Boland)  [to  appear,  SIAM  J.  Num.  Analysis]. 

3.  Solution  technique  for  discrete  incompressible  Navier-Stokes 
equations  [in  preparation]. 

4.  Gaussian  elimination  with  noninvertible  pivots  (with  M.  D. 
Gunsburger)  [to  appear,  J.  Linear  Algebra]. 


